Bounds in Probability

Adapted from: SOSTOOLS' SOSDEMO8 (See Section 4.8 of SOSTOOLS User's Manual)

The probability adds up to one.

μ0 = 1

The mean is one.

μ1  = 1

The standard deviation is 1/2.

σ = 1/2

0.5

The second moment E(x^2) is:

μ2 = σ^2 + μ1^2

1.25

We define the moments as follows:

using DynamicPolynomials
@polyvar x
monos = [1, x, x^2]
using SumOfSquares
μ = measure([μ0, μ1, μ2], monos)

MultivariateMoments.Measure{Float64, DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}([1.25, 1.0, 1.0], DynamicPolynomials.Monomial{true}[x², x, 1])

We need to pick an SDP solver, see here for a list of the available choices. We use SOSModel instead of Model to be able to use the >= syntax for Sum-of-Squares constraints.

using CSDP
solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
model = SOSModel(solver);

We create a polynomial with the monomials in monos and JuMP decision variables as coefficients as follows:

@variable(model, poly, Poly(monos))

\[ (_[1])x^{2} + (_[2])x + (_[3]) \]

Nonnegative on the support:

K = @set 0 <= x && x <= 5
con_ref = @constraint(model, poly >= 0, domain = K)

\[ (_[1])x^{2} + (_[2])x + (_[3]) \text{ is SOS} \]

Greater than one on the event:

@constraint(model, poly >= 1, domain = (@set 4 <= x && x <= 5))

\[ (_[1])x^{2} + (_[2])x + (_[3] - 1) \text{ is SOS} \]

The bound (we use LinearAlgebra for the ⋅ syntax for the scalar product):

using LinearAlgebra
@objective(model, Min, poly ⋅ μ)

\[ 1.25 {\_}_{1} + {\_}_{2} + {\_}_{3} \]

We verify that we found a feasible solution:

optimize!(model)
primal_status(model)

FEASIBLE_POINT::ResultStatusCode = 1

The objective value is 1/37:

objective_value(model)

0.027027027670456782

The solution is (12x-11)^2 / 37^2:

value(poly) * 37^2

\[ 144.00007395188794x^{2} - 264.000359545551x + 121.00026798654635 \]

This page was generated using Literate.jl.